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-In  many  surface  water  bodies,  water  temperature  closely  follows  ambient  air 
temperature.  This  means  that  warmer  water  in  winter  absorbs  hast  from  below. 
The  extent  and  pattern  of  winter  heat  gain  la  constrained  by  the  fact  that 
tha  water  temperature  does  not  fall  below  the  f reeling  point.  On  the  basis  of 
a  few  simple  assumptions,  governing  equations  are  solved  hers  pertaining  to 
heat  flow  In  bottom  sediments.  The  rneults  are  presented  in  general  nondi- 
manalonal lted  curves.  These  allow  estimation  of  water/eedlment  hast  flux  — ^ 
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(or  any  particular  caae,  given  truncation  of  tha  water  taaiparature  curve  at  the 
f reeling  point.  The  uaer  muat  aupply  pertinent  yearly  air  temperature  mean 
and  amplitude  of  variation,  together  with  the  thermal  dlffualvity  for  the  bottom 
material.  The  governing  equationa  are  eolved  ualng  a  higher  order  finite 
element  method  which  eolvea  directly  for  temperature  gradlenta  and  hence  for 
heat  flux.  Thue  the  method  provldea  particularly  accurate  flux  valuee  at  high 
efficiency.  The  reeulta  llluatrate  in  detail  how  winter  water  heat  gain  la  leaa 
in  caaea  where  mean  air  temperaturea  are  lower. 
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This  report  mss  prspsrsd  by  Dr.  Ksvin  O'Neill  of  tha  Gsotschnlcsl 
Rssssreh  Branch,  Experimental  Engineering  Division,  and  Dr.  Gaorge  Ashton 
of  tha  Snow  and  lea  Branch,  laaaarch  Division,  U.S.  Army  Cold  Regions  Re¬ 
search  and  Engineering  Laboratory.  Tha  study  was  supported  by  tha  Office, 
Chief  of  Engineers,  Directorate  of  Civil  Works,  under  Work  Unit  31594, 
Formation  and  Breakup  of  Reservoir  Ice  Cover  and  Effects  on  Thermal  Energy 
Budaet  Computations. This  work  unit  is  part  of  the  Environmental  and  Water 
Quality  Operational  Studies  program  managed  by  Jerome  L.  Mahloch  of  the 
U.S.  Army  Engineer  Waterways  Experiment  Station. 

The  report  was  technically  reviewed  by  Dr.  Dennis  Ford  of  WES. 

The  contents  of  this  report  are  not  to  be  used  for  advertising  or  pro¬ 
motional  purposes.  Citation  of  brand  names  does  not  constitute  an  official 
endorsement  or  approval  of  the  use  of  such  commercial  products. 
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INTRODUCTION 


One  of  the  components  of  the  total  energy  budget  of  a  water  body 
la  the  energy  stored  in  and  released  from  the  bottom  materials  (usually 
sediments).  The  heat  flux  is  seasonal,  with  heat  conducted  into  the  bottom 
during  summer  when  the  water  is  warm  and  released  to  the  cold  water  during 
winter.  This  paper  explores  the  magnitude  of  this  heat  flux,  particularly 
in  winter,  when  the  water  temperature  Is  truncated  at  0°C  from  the  approxi¬ 
mately  sinusoidal  variation  that  would  otherwise  persist. 

NATURE  OF  THE  PROBLEM 

A  reasonable  representation  of  the  annual  temperature  of  a  water 
body  is  a  sinusoidal  variation  more  or  less  following  the  average  air  tem¬ 
peratures.  However,  when  the  air  temperature  falls  below  0*C  the  water 
temperature  cannot  follow  it,  because  of  the  state  condition;  heat  loss 
from  water  at  0°C  to  the  subfreezing  atmosphere  manifests  itself  in  ice 
formation  while  the  liquid  water  remains  at  0*C.  In  lakes  with  a  complete 
ice  cover,  the  cover  also  seals  the  water  body  from  energy  exchange  with 
the  atmosphere.  Changes  in  the  water  temperature  are  dominated  by  heat 
flux  from  the  bottom  materials,  which  gradually  raises  the  water  tempera¬ 
ture  through  the  winter.  While  bottom  heat  fluxes  in  summer  are  generally 
not  important  components  of  the  energy  budget  of  water  bodies,  they  are 
important  during  the  winter. 

A  completely  rigorous  treatment  of  the  bottom  heat  flux  contribution 
requires  specification  of  the  water  temperature  at  the  sediment/water 
Interface  to  include  the  gradual  warming  due  to  the  heat  flux.  Neverthe¬ 
less,  a  good  approximation  is  to  assume  a  constant  bottom  surface  tempera¬ 
ture  of  0"C  during  subfreezing  weather  (Fig.  I).  We  also  note  that  the 
constant  0®C  water  temperature  during  the  winter  is  an  excellent  approxima¬ 
tion  for  rivers  and,  while  the  associated  heat  flux  causes  only  a  small 
temperature  rise  (on  the  order  of  0.01*C),  the  heat  flux  itself  may  have  a 
significant  effect  on  frazil  ice  in  the  flow.  Finally,  we  note  that  the 
strictly  sinusoidal  surface  temperature  variation  on  a  semi-infinite  medium 
is  a  classic  problem  in  conductive  heat  transfer  and  has  been  solved 
analytically 


Figure  1.  Generalized  river  bottom 
temperature  variation  with  time. 
Between  time  t|  and  t2  the  ambient 
air  temperature  follows  the  dashed 
line,  while  the  water  temperature  does 
not  drop  below  0*C. 


1 Cars law,  H.S.  and  J.C.  Jaeger  (1973)  Conduction  of  Heat  in  Solids.  Second 
Edition.  London:  Oxford  University  Press,  p.  64-70. 
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ANALYSIS 


The  specified  surface  temperature  variation  Is  shown  in  Figure 
1.  Tg  is  the  sediment  surface  temperature  at  any  time  t,  Tg  is  the 
mean  and  Ta  the  amplitude  of  the  strictly  sinusoidal  (air)  surface 
temperature;  T  is  the  mean  of  the  truncated  sinusoidal  variation  and  is  the 
temperature  which  will  be  approached  at  very  great  depths  in  the  bottom 
after  imposition  of  the  truncated  surface  temperature  over  many  years. 

These  quantities  ideally  would  be  obtained  from  several  years  of  water 
temperature  measurements.  As  a  practical  matter,  however,  a  reasonable 
approximation  is  to  fit  a  sinusoidal  variation  to  average  monthly  air 
temperatures  obtained  from  nearby  meteorological  stations,  since  the  water 
temperature  variations  of  most  shallow  water  bodies  closely  follow  air 
temperature  variations. 


The  governing  heat  conduction  equation  is 


II 

3t 


(1) 


where  T(x,t)  is  the  temperature  in  the  bottom  at  distance  x  below  the 
sediment/water  interface  and  K  ■  k/pc  where  k  is  the  thermal  conductivity 
of  the  bottom  material,  p  is  the  density,  and  c  is  the  heat  capacity.  (The 
Celsius  temperature  scale  is  assumed  below.  In  principle,  any  other  scale 
would  result  in  identical  expressions,  as  long  as  the  freezing  temperature 
of  water  is  taken  as  zero  degrees.) 


The  boundary  conditions  are 

T(0,t)  *  1  +  T  sin  2*ft,  t  not  between  t.  and  t.  (2) 

m  a  i  £ 


T(0,t)  -  0  *1  <  c  <  *2  (3) 

T(-,t)  -  T  (4) 


T(x,0)  -  T.  (5) 


In  eq  2,  f  is  the  frequency,  which  we  take  to  be  1  year"*. 

It  is  convenient  to  non-dimens lonallze  the  equation  by  introducing 

t-T 

e  ■  J  ■“  so  at  x  -  0,  1  <  8  <_  -  r  (6) 

a  8 

where  r  ■  Ta/Ta.  Also,  using 

t  ■  ft,  $  ■  x//K/f  (7) 

the  equations  assusw  the  form 
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T  -  T 

8(»,t)  "  8(5.0)  -  -J-S.  (9) 

a 

and  the  boundary  condition  in  eq  2  and  3  is  scaled  in  the  same  manner. 

Since  we  are  primarily  Interested  in  the  temperature  gradient  at  the 
sediment/water  interface,  we  seek  to  determine  39/35  there  as  a  function  of 
r  and  time.  At  this  point  the  advantages  of  nondioensionallzatlon  are 
apparent:  all  possible  results  may  be  displayed  in  a  single  set  of  curves 
in  30/35  versus  t«  Once  38/35  is  determined,  the  temperature  gradient  at 
the  sediment /water  interface  is  determined  by  the  transformation 

JT  _  ii  -  t  3C  36  _  Ta  30  nov 

3x  a  3x  a  3x  35  35  * 

The  problem  was  analyzed  using  a  finite  element  formulation  described 
in  more  detail  in  the  next  section.  The  resulting  output  for  the  case  of 
r  -  1  (no  truncation)  was  compared  to  the  analytical  solution  and  matched 
to  three  significant  figures.  When  r  *  1,  the  program  was  run  for  about 
five  periods  to  remove  transients.  About  20-36  time  steps  per  period 
(1  year)  were  used  with  a  10-element  mesh  (11  mesh  points),  which  resulted 
in  a  very  low  cost  per  run.  This  efficiency  was  due  to  1)  a  high  order  of 
interpolation  between  mesh  points  and  2)  the  particular  Implicit 
formulation  in  time. 


THE  FINITE  ELEMENT  ANALYSIS 

In  this  section,  a  brief  outline  of  the  basics  of  the  Galerkln 
finite  element  method  is  given.  It  is  by  no  means  sufficient  for 
answering  all  questions  concerning  the  method  and  ’'he  reader  is  referred  to 
the  voluminous  literature  for  further  information2. 

In  the  finite  element  method,  mathematical  functions  are  expressed 
as  sums  of  finite  series,  utilizing  a  preselected  set  of  so-called 
basis  functions.  Thus,  if  F  is  the  function  in  question,  it  is 
expressed  as 


F 


N 

-  I  Fj  U< 

J-l  3  3 


(ID 


2E.g. :  Zienkiewlcz,  O.C.  (1977)  The  Finite  Element  Method.  Third  Edition. 
New  York:  McGraw-Hill. 
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where  the  F.  are  coefficients,  N  is  a  finite  integer,  and  the  functions 
Uj(x)  are  basis  functions.  In  heat  flow  problems,  the  coefficients  are 
usually  time-dependent,  varying  according  to  the  particulars  of  the 
problem,  whereas  the  basis  functions  are  written  in  space  in  a  temporally 
invariant  manner: 


N 

F  -  l  F  (t)  U  (x)  .  (12) 

j-1  J  J 

The  basis  functions  may  be  of  many  kinds,  but  in  any  case,  are  con¬ 
ceptually  the  first  N  members  of  a  set  which,  in  its  entirety,  could  theo¬ 
retically  express  the  function  F  exactly.  The  use  of  such  a  series  and  the 
subsequent  manipulations  of  it  are  reminiscent  of  the  use  of  Fourier  series 
in  solving  many  problems.  A  major  difference  is  that  the  sinusoidal  func¬ 
tions  in  a  Fourier  series  are  in  general  nonzero  throughout  the  entire 
domain.  In  contrast,  basis  functions  for  finite  elements  are  typically 
locally  defined,  each  being  nonzero  only  within  the  space  elements  adjacent 
to  some  particular  mesh  point  (‘‘node")  with  which  the  basis  function  is 
associated  (Fig.  2).  An  arbitrary  function  may  be  approximated  using  the 
functions  depicted  in  Figure  2  by  selecting  appropriate  coefficients,  in 
this  case  the  approximated  function's  values  at  the  node  points  (Fig.  3). 

The  basis  functions  need  not  be  linear,  nor  must  the  coefficients  cor¬ 
respond  to  the  function's  values  only,  but  may  correspond,  for  example,  to 
its  derivatives. 

Consider  now  a  differential  equation  to  be  solved,  written  in 
the  abstract  as 


UyJ  -  0 


Figure  2.  Typical  linear  finite 
element  basis  functions,  on  a  uniform 
grid.  U5,  for  example,  rises  to  a 
magnitude  of  1  at  node  point  X5, 
declines  to  0  at  X4  and  xg,  and  is 
uniformly  0  everywhere  except  between 
X4  and  xg.  The  domain  between  any 
two  node  points  constitutes  one  finite 
element. 


(13) 


Figure  3.  A  function  F,  shown  as  the 
solid  line,  may  be  approximated  in  the 
manner  of  eq  12.  Using  the  basis 
functions  in  Figure  2  with  the  function 
values  of  F  at  the  nodes  as  the 
coefficients,  one  obtains  the 
approximation  shown  as  the  dashed 
line. 
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where  L  Is  some  differential  operator.  For  example,  if  L  corresponds  to 


then  eq  13  denotes 


L  .  i_  .  i_ 
3T  C2  ’ 


iZ  _ 

3T  a*2 


The  solution  to  the  equation  is  some  function  y,  for  which  the  rela¬ 
tion  expressed  in  eq  13  holds  exactly.  If,  instead  of  y,  one  substitutes 
into  eq  13  an  approximation  of  it,  denoted  y,  then  the  equality  will  not  be 
exactly  satlslfed,  and 

L[y]  -  R(x)  .  (16) 


R  is  called  the  "residual,”  and  inasmuch  as  y  does  not  satisfy  eq  13 
exactly,  R  will  not  be  uniformly  zero. 

It  is  the  undertaking  of  a  class  of  numerical  methods  called 
"weighted  residual  methods”  to  find  approximations  y  which  come  desirably 
close  to  satisfying  the  governing  equation.  Specifically,  it  is  required 
that 

/  w  R  dx  «  J  w  L[y]  dx  -  0  (17] 

S  S 

where  w^(x)  is  one  of  a  selected  set  of  weighting  functions,  and  S  is 
the  domain  considered.  For  example,  if  w^  is  uniformly  equal  to  one, 
then  eq  17  requires  only  that  the  net  or  sum  of  the  residual  over  the 
entire  domain  be  zero  (an  unreliable  criterion).  The  functions  w^  may  be 
any  of  a  considerable  variety  of  possibilities,  each  giving  rise  to  a 
different  specific  method  of  attack,  such  as  the  collocation  method,  the 
subdomain  method,  and  the  least  squares  method. 

In  the  Galerkln  finite  element  method,  the  set  of  weighting  functions 
is  chosen  to  be  the  basis  functions  themselves.  Both  experience  and 
theoretical  analysis  justify  the  expectation  that  the  use  of  the  as 
weighting  functions  can  produce  accurate  approximations  to  the  desired  y. 
Thus,  one  approximates  y  as 


y  -  l  Y  (t)  0  (x)  ( 

j-1  3  3 

where  the  Yj  are  initially  known,  and  in  accordance  with  eq  17  one  then 
writes 


/  \Ji  U^YjUj]  dx  -  0 
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or 


N 

l  {/  «  L(U  ]dx}  Y  -  0.  (20) 

j-1  S  J  J 

One  such  equation  results  for  each  in  the  set  of  N.  In  matrix 
notation,  eq  20  is  equivalent  to 

[A] {y}  -  0.  (21) 

When  the  boundary  conditions  are  incorporated,  eq  21  constitutes  a  set  of  N 
linear  algebraic  equations  in  the  N  coefficients  sought,  which  can  be 
solved  by  any  appropriate  standard  method. 

All  results  reported  below  were  obtained  using  a  particular  set 
of  basis  functions  known  as  Hermite  basis  functions.  These  functions 
have  the  property  of  generating  a  particularly  accurate,  higher  order 
(cubic)  interpolation  of  functions  between  the  node  points.  They  are 
also  selected  so  that  half  the  unknown  coefficients  correspond  to 
function  values  at  the  nodes;  the  other  coefficients  correspond  to  the 
function's  gradient  values.  Thus  gradients  of  the  unknown  function  at 
the  node  points  are  obtained  directly  in  the  course  of  solving  the 
problem,  without  further  manipulation  or  numerical  differentiation. 

This  system  is  especially  desirable  for  the  problem  treated  here, 
since  flux  values  based  on  solution  gradients  are  the  result  ultimately 
sought. 

RESULTS  AND  INTERPRETATION 

In  Figure  4  are  presented  curves  of  36/35  at  x  -  0,  for  r  -  1.0  (no 
truncation),  0.9,  0.8,  0.6,  0.4  and  0.0  (mean  air  temperature  *  0#C).  It 
is  clear  from  Figure  4  that  the  truncation  of  the  surface  temperature  has  a 
significant  effect  on  the  temperature  gradient,  particularly  during  the 
winter. 

Practical  use  of  these  results  to  calculate  the  bottom  heat  flux 
contribution  requires  knowledge  of  the  surface  temperature  mean  and 
amplitude,  and  the  bottom  thermal  diffusivity  K.  For  shallow  water  bodies 
it  is  reasonable  to  use  mean  monthly  air  temperature  to  determine  the 
former.  The  thermal  diffusivity  depends  upon  the  thermal  properties  of  the 
bottom  materials.  However,  the  bottom  materials  of  most  water  bodies  are 
sediments  saturated  with  water.  Measurements  of  the  thermal  conductivity, 
from  which  the  diffusivity  may  be  calculated,  show  that  bottom  materials 
may  be  considered  to  have  the  same  diffusivity  as  still  water3. 

Accordingly,  in  the  examples  which  follow,  the  diffusivity  K  will  be  taken 
as  1.31  x  10“  ^  m^s~l,  which  is  the  value  for  water  at  0#C. 

As  an  example,  we  will  use  the  mean  air  temperatures  for  Hanover,  New 
Hampshire  (Fig.  5).  The  mean  annual  temperature  (Ta)  is  6.3*C  and  the 


%cGaw ,  R.  (1974)  Thermal  conductivity  of  organic  sediments  from  two 
Wisconsin  lakes.  U.S.  Army  Cold  Regions  Research  and  Engineering 
Laboratory  Special  Report  129. 


»  : 


6 


Figure  4.  Nondlmena tonal  river  bottom 
temperature  gradient,  evaluated  at  the 
sedlment/water  Interface,  as  a  function 
of  time  and  r.  The  nondlmensional 
river  bottom  temperature  8S  is  also 
shown.  In  the  example  in  the  text,  r  - 
0.444,  and  the  8S  curve  is  truncated 
at  the  level  shown. 


Figure  S.  An  example  using  the  curves 
in  Figure  4  for  Hanover,  N.H. 

Variation  of  monthly  mean  temperatures 
is  shown  at  the  top,  based  on  data  from 
nearby  weather  stations.  The 
corresponding  river  bottom  heat  flux 
response  is  shown  below. 


amplitude  (Ta)  is  14.2°C,  from  which  r  -  6.3/14.2  -  0.444.  Initial  time 
for  the  sinusoid  that  fits  the  air  temperature  variation  is  approximately 
mid-April.  Corresponding  heat  flux  response  over  the  months  may  then  be 
read  off  Figure  4,  with  values  between  the  curves  for  r  *  0.4  and  r  •  0.6. 
The  transformation  to  dT/3x  from  36/3£  i« 


JT  m  Ta  38  m  _ 14.2*C _  j)8 

**  /it7f  35  (I.31xl0_7m2s”1  .  86400  .  365  s)1/2  3^ 
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The  actual  heat  flux  Is  then  found  from  It  aT/jx  where  k  is  the  conduc¬ 
tivity  of  the  bottom  material  (taken  as  the  value  for  water  at  0®C;  k  ■ 

0.55  W  Figure  5  shows  heat  flux  values  based  on  this  value 

of  k  and  on  values  from  Figure  4.  Examination  of  Figure  5  shows  that  the 
heat  gain  to  the  water  reaches  a  maximum  in  November  (about  the  time  of 
initial  ice  cover  formation)  and  then  decreases,  rapidly  at  first,  over  the 
winter. 

This  analysis  is  most  applicable  to  rivers  and  to  reservoirs  where 
sufficient  through-flow  exists  during  winter  to  ensure  that  the  bottom 
water  is  at  0°C  throughout  the  winter.  In  natural  lakes  with  little 
through-flow,  the  Initial  temperature  at  the  onset  of  ice  is  often  somewhat 
above  0°C  (but  seldom  above  4®C).  In  this  case  a  good  first  approximation 
is  to  use  the  same  procedure  but  with  an  r  value  determined  using  the 
initial  temperature,  rather  than  0®C.  However,  since  the  initial  tempera¬ 
ture  at  onset  of  ice  is  determined  by  the  particular  sequence  of  meteoro¬ 
logical  events  after  the  4°C  isothermal  state  has  been  achieved,  and  these 
events  vary  from  year  to  year,  it  is  difficult  to  predict  the  initial  temp¬ 
erature. 

As  was  discussed  earlier,  a  complete  solution  to  the  problem  would 
include  the  effect  of  elevation  of  the  bottom  water  temperatures  as  a 
result  of  the  heat  flux  from  the  bottom.  The  effect  of  this  warming  would 
be  to  cause  the  falling  limb  of  the  “winter  heat  gain“  curve  (Fig.  5)  to 
fall  further,  but  would  have  little  effect  on  the  higher  values  of  heat 
flux  near  the  beginning  of  the  winter. 
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